using System;
using System.Reflection;
using NBody;
using NBody.Analysis;
using NBody.Cosmology;
using OptionParser;
using DataFile;

[assembly: AssemblyTitle ("NBodyAngle")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Calculate the change in azimuthal angle as a function of radius.")]
[assembly: AssemblyCopyright ("2006 Joseph D. MacMillan")]

public class AngleOptions : InputOutputGetOptions
{
    [Option ("-b", "Number of radial bins")]
    public int NumBins = 100;
    
    [Option ("-e", "Softening length")]
    public double Softening;
    
    [Option ("-f", "File for potential data")]
    public string PotFile;
}

public class NBodyAngle
{
    public static void Main(string[] args)
    {
        // Process command-line options
        AngleOptions opts = new AngleOptions();
        opts.ProcessArgs(args, "angle");
                        
        // Read in a N-body system from the given file
        NBodySystem s = NBodySystem.Read(opts.InFile);
            
        if (opts.PotFile != "")
            Energy.SetPotentialFromData(ref s, Table.Read(opts.PotFile));
        else
            Energy.CalculatePotentialShell(ref s, opts.Softening);
        
        double e_tot;
        Energy.CalculateTotalEnergy(ref s, out e_tot);   
        
        // now for each shell, calculate delta phi
        s.Sort();
        
        Console.Error.WriteLine("Binning in radial shells with logarithmic spacing");
        double rmax = s.Rmax + 1e-10;
        double rmin = s.Rmin;
            
        double deltaR = Math.Log10(rmax / rmin) / (double)opts.NumBins;
        
        Table phi = new Table(opts.NumBins, 4, "Radius, rms E, rms j, delta phi");
        phi.Clear();
        
        // first calculate rms E and j values
        int[] num = new int[opts.NumBins];
        foreach (EnergyParticle p in s)
        {
            int pos = (int)((Math.Log10(p.Radius) - Math.Log10(rmin)) / deltaR);
            num[pos]++;
            phi[pos, 1] += p.Energy;
            phi[pos, 2] += Math.Sqrt(p.J2);
        }
             
        for (int i = 0; i < opts.NumBins; i++)
        {
            phi[i, 0] = Math.Pow(10.0, Math.Log10(rmin) + 0.5 * deltaR + deltaR * i);
            if (num[i] > 0)
            {
                phi[i, 1] /= (double)num[i];
                phi[i, 2] = (phi[i, 2]) / (double)num[i];
            }
        }
                        
        double[] r = new double[s.NumParts];            
        double[] g1 = new double[s.NumParts];         
        for (int i = 0; i < s.NumParts; i++)
            r[i] = s[i].Radius;
            
        for (int i = 0; i < opts.NumBins; i++)
        {
            double e = phi[i, 1];
            double jrms = phi[i, 2];
            
            for (int j = 0; j < s.NumParts; j++)
            {
                g1[j] = 2.0 * (e - ((EnergyParticle)s[j]).Potential) - jrms*jrms / r[j] / r[j];
                if (g1[j] < 0.0)
                    g1[j] = 0.0;
                else
                    g1[j] = 1.0 / Math.Sqrt(g1[j]) / r[j] / r[j];
            }
            phi[i, 3] = 2 * jrms * Integration.Trapezoidal(r, g1, 0, s.NumParts - 1);
            phi[i, 0] = Math.Log10(phi[i, 0]);
        }
        
        phi.Print(opts.OutFile);            
        
    }
}
